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ABSTRACT 

' We describe an implicit general relativistic hydrodynamics code. The evolution 

equations are formulated in comoving coordinates. A conservative finite differencing 
^ • of the Einstein equations is outlined and artificial viscosity and numerical diffusion are 

<^ . discussed. The time integration is performed with AGILE, an implicit solver for stiff 

^ , algebro-differential equations on a dynamical adaptive grid. We extend the adaptive 

CN [ grid technique, known from nonrelativistic hydrodynamics, to the general relativistic 

I application and identify it with the concept of shift vectors in a 3+1 decomposition. The 

K*" ' adaptive grid minimizes the number of required computational zones without compro- 

■ mising the resolution in physically important regions. Thus, the computational effort 

lO ' is greatly reduced when the zones are subject to computationally expensive additional 

processes, such as Boltzmann radiation transport or a nuclear reaction network. We 
present accurate results in the standard tests for supernova simulations: Sedov's point 
blast explosion, the nonrelativistic and relativistic shock tube, the Oppenheimer-Snyder 
dust collapse, and homologous collapse. 
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^ ' Subject headings: supernovae: general — hydrodynamics — gravitation — relativity — shock 

^ ■ waves — methods: numerical 



1. Introduction 

Consider a physical variable y that evolves in time t according to a differential equation dy /dt = 
f{y). The simplest implicit finite differencing relates the value of at time to its previous value 
y""-*^ at time f^'^ by 

n _ n-1 

- — = fiy"")- (1) 

fn _ ^n— 1 J \f J \ / 

This equation can be solved for at each time step n = 1 . . . oo. Depending on the properties of 
the function /, a computationally expensive solution of a nonlinear algebraic equation is required 
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in most cases. This is not the case with the expUcit finite differencing where the evaluation of / 
is based on the old values y"^~^. The two discretization variants appear to be similar for small 
time steps where — is small. However, the behavior dramatically differs for larger 

steps: If a physical equilibrium with the static solution = f{y°°) exists, one easily checks with 
equation (1) that the equilibrium solution can be approached in one large time step because of 
the bound value of \y°° — Explicit finite differencing, on the other hand, does not lead 

to a physically meaningful solution with time steps larger than a characteristic time scale that is 
determined by the function /. Implicit methods guarantee smooth transitions between dynamical 
phases and equilibrium phases. On a long time scale, they automatically adjust the equilibrium to 
exterior dependencies and provide an averaging approximation where physical oscillations around 
the equilibrium are not of primary interest. Although implicit schemes certainly are capable of 
tracking the full dynamical solution on short time scales, a calculation with time steps smaller than 
the characteristic time scale is much more efficiently done with explicit finite differencing. The 
reduction in the number of time steps with the implicit technique has to compensate the additional 
effort invested with the solution of the algebraic equations. 

Processes in astrophysics involve characteristic time scales that differ by orders of magnitude. 
For example Henyey ct al. (1959) had to introduce implicit finite differencing in stellar evolution 
calculations, where nuclear reactions as well as hydrodynamics were assumed to be quasi-static on 
the time scale of stellar evolution. The system of algebraic equations is traditionally solved with 
a Newton-Raphson scheme. The numerical effort increases steeply with the number of unknowns 
that have to be determined in one time step. Thus, implicit approaches are predestined for systems 
of moderate size where the evolution time scale is comparable to the time scale of the underly- 
ing physical processes somewhere in space or time, and exceeds the latter by far elsewhere or at 
another time. Because of this size restriction, implicit hydrodynamics has mainly been applied 
in one-dimensional systems, such as in spherically symmetric hydrodynamics codes for supernova 
simulations (Schinder & Bludman 1989; Swesty 1995; Yamada 1997). Spherical symmetry addi- 
tionally allows the use of Lagrangian comoving coordinates and a full treatment of general relativity 
(Misner &: Sharp 1964). The evolution of the metric can be included self-consistently into the so- 
lution for the nonlinear algebraic equations that emerge from the implicit finite differencing. 

Fig. (1) shows the connection between several technical features that go along with a hydro- 
dynamics code that is based on implicit finite differencing. Eggleton (1971) used an adaptive grid 
in his quasistatic stellar evolution code. Together with the hydrostatic equations, a grid equation 
is solved that determines where a fixed number of grid points is located to optimally sample the 
profile of the star. A dynamical adaptive grid was suggested by (Harten k, Hyman 1983) and the 
formulation of radiation hydrodynamics on a dynamical adaptive grid was provided by Winkler, 
Norman, & Mihalas (1984) and implemented in Newtonian (Winkler & Norman 1986) and special 
relativistic (Norman & Winkler 1986) hydrodynamics. Dorfi & Drury (1987) found a simple and 
robust implicit prescription to guide the dynamics of the adaptive grid according to the evolution of 
the physical state. Note the difference between the dynamical adaptive grid where a fixed number of 
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grid points continously moves through the computational domain and a static adaptive grid where 
static grid points are inserted or removed according to resohition requirements. The dynamical 
adaptive grid technique found applications in astrophysical high-resolution calculations (Dorfi & 
Gautschi 1989; Mair 1990), and a well documented implementation of the dynamical adaptive 
grid technique for Newtonian hydrodynamics with radiative transfer in different one-dimensional 
geometries has been provided by Gehmeyr &; Mihalas (1994) in the code TITAN. 

We resume the description of Fig. (1) with the relationship between a dynamical adaptive 
grid and implicit hydrodynamics. The solver of the algebraic equations easily incorporates the 
solution of the implicit grid equation. The zones are dynamically concentrated in the regions 
of action or where a high resolution is required. This assures that most of the zones actively 
contribute to the interesting part of the evolution and do not clump together in regions of well- 
known physical equilibrium. The dynamical adaptive grid can reduce the required number of grid 
points by orders of magnitude with respect to an equidistant grid, and the implicit scheme greatly 
benefits from this size reduction in the solution vector. Since the dynamical grid changes neither 
the total number of grid points nor its ordering, the structure of the Jacobian that has to be 
inverted in the Newton-Raphson scheme remains unchanged during the simulation. Furthermore, 
the high resolution achieved with the adaptive grid is only advantageous in combination with 
implicit hydrodynamics because the Courant-Friedrichs-Lewy condition (Courant, Friedrichs, & 
Lewy 1928) would impose severe restrictions on the maximum time step in explicit schemes. 

Of course, the adaptive grid also comes with difficulties. The moving grid relates different 
spatial points in the physical state vector across one and the same time step. This requires a 
proper handling of advection in the formulation of the physical equations on the adaptive grid. 
Moreover, the advection introduces extended coupling between neighboring zones and therewith 
increases the density of the Jacobian in the Newton-Raphson scheme. A conservative formulation 
of the physical equation provides a solid basis for the discretization of the advection. Conserved 
quantities are transfered from one zone into a neighboring zone. A formulation of the dynamics 
in terms of conservation laws has additional benefits: (i) Fundamental conservation laws are also 
valid at discontinuities and allow an accurate numerical solution - e.g., for the propagation of shock 
waves. This is not the case with arbitrary finite difference approximations, (ii) The integration of 
a conserved quantity over adjacent fluid elements does not depend on the fluxes at the enclosed 
boundaries. Thus, discretization errors have mainly local influence. This advantage is important 
for the implicit solution of problems involving scale differences of many orders of magnitude, (iii) 
The integration over the whole computational domain of the single-fluid-element conservation laws 
leads to the total conservation of fundamental physical quantities, independent of the resolution in 
the conservative flnite difference representation. These contributions of a conservative formulation 
to the adaptive grid and implicit solution are also drawn in Fig. (1). 

In this paper, we describe our implementation of the dynamical adaptive grid technique in 
general relativistic hydrodynamics. The dynamical adaptive grid is nothing other than a continuous 
and time dependent shift in the coordinates with respect to Lagrangian or Eulerian coordinates. 
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Fig. 1. — The interdependence between general relativity, implicit finite differencing, conservative 
equations, and the dynamical adaptive grid in spherically symmetric hydrodynamics. 
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This freedom of coordinate choice has long been explored in the general relativistic case, where 
in addition to the time slices, shift vectors have to be defined in a 3+1 decomposition (Arnowitt, 
Deser, & Misner 1962; Smarr &: York 1978). There is no difference between a dynamical adaptive 
grid and a specific choice of shift vectors in the general relativistic description. However, as far 
as the numerics is concerned, the formulation with shift vectors only provides a set of analytical 
equations, while the Newtonian adaptive grid technique can be generalized to a detailed prescription 
for a stable numerical implementation of shift vectors in the general relativistic case. Thus, instead 
of directly applying a discretization of appropriate shift vectors, we first search for a conservative 
formulation of hydrodynamics in Lagrangian coordinates and then apply the generalized adaptive 
grid technique to implement the shift vectors that lead to an optimal resolution of the physical 
state at important locations. 

All the features drawn in Fig. (1) come together in an important astrophysical application: 
core collapse supernovae. The compact object at the center of the star requires a general relativistic 
treatment (Bruenn, DeNisco, &; Mezzacappa 2001; Liebendorfer et al. 2001). The time scales of 
radiation transport, hydrodynamics, and neutrino diffusion differ by orders of magnitude and call for 
implicit techniques. The computational domain covers densities from 10^ g cm~^ in the outer layers 
of the calculational domain to several times 10^^ g cm^"^ in central regions of the proto-neutron star. 
At these higher-than-nuclear densities the stiffness of the equation of state is responsible for very 
short hydrodynamical time scales. The speed of sound waves reaches several 10% of the velocity 
of light. If one is not interested in the ringing of the neutron star, this time scale is prohibitive 
for the calculation of the long term physics in the outer range with schemes implementing explicit 
finite differencing. The supernova shock formed after stellar core collapse and bounce requires 
reasonable resolution. Moreover, in the case of a supernova explosion, the adaptive grid allows a 
smooth propagation of the shock through the outer layers; and in the case of black hole formation, 
a robust implementation of shift vectors might help to shift the arising coordinate singularity at the 
Schwarzschild horizon out of the computational domain (Liebendorfer &; Thielemann 2000). Even 
the conservative formulation of the hydrodynamic equations is physically motivated: A meticulous 
energy conservation in the finite difference scheme is required to allow statements on the success 
of a supernova explosion (Mezzacappa et al. 2001). The balance of gravitational, internal, and 
neutrino energy is two orders of magnitude smaller than the individual constituents. The total 
energy is then comparable to the explosion energy, i.e. about one percent of the energy stored in 
the neutrino field. The bookkeeping of this energy has to be based on a careful implementation of 
the mutual energy and momentum transfer between the neutrino radiation field and the matter. For 
the determination of the explosion energy to an accuracy of at least one order of magnitude, one may 
tolerate an energy drift of not exceeding a fraction of 10~^ of the internal, gravitational or neutrino 
radiation energy. This drift extends over the whole calculation time of more than 10000 time steps, 
so that systematic deviations of energy conservation may not exceed a fraction of 10~^ per time 
step. On current computer hardware, it is not yet possible to simulate in multi-dimensional general 
relativity the confluence of accurate multigroup neutrino transport, high resolution hydrodynamics, 
magnetic fields, and up-to-date nuclear physics input. The complexity of supernova explosions is 
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best explored in separate restricted investigations. Spherically symmetric simulations are able to 

incorporate the best available neutrino transport techniques (Rampp & Janka 2000; Mezzacappa 
et al. 2001) in self-consistent simulations with high resolution and allow the detailed study of high 
density weak interactions and nuclear physics. 

In the following section, we write down conservative equations of general relativistic hydro- 
dynamics in spherical symmetry and discuss the relation between shift vectors and the adaptive 
grid. The adaptive grid technique of Winkler, Norman, k, Mihalas (1984) is extended to general 
relativity and applied in combination with the grid equation of Dorfi k, Drury (1987). We then 
proceed with a description of our approach to automatically and implicitly integrate a system of 
stiff ordinary differential equations in time and outline the spatially staggered finite differencing of 
the physical equations. Numerical diffusion and artificial viscosity are discussed. In section 4 we 
run the standard supernova test simulations with our code, AGILE, and compare the results to 
analytical solutions and the Piecewise Parabolic Method (PPM) (Colella &; Woodward 1984). 

2. Conservative Equations on the Dynamical Adaptive Grid 

2.1. Relativistic hydrodynamics in spherical symmetry 

General relativistic four-dimensional space-time can be described in terms of consecutive three- 
dimensional time slices (Arnowitt, Descr, &: Misner 1962). The definition of the time slices and 
the choice of coordinates within the slices is not given by the dynamics of the matter; it has to 
be imposed from the outside as a part of the description. The time slices comprise the events 
that are made simultaneous with respect to coordinate time. They are defined by a lapse function 
that describes the proper time lapse to the infinitcsimally close neighbor time slice t + 8t along a 
path orthogonal to the time slice at time t. The three-dimensional coordinates within a time slice 
are then specified by a shift vector function describing the spatial shift of points with the same 
coordinate label between infinitcsimally close neighbor slices. Reference for the shift is again the 
path orthogonal to the shce (Smarr & York 1978). 

In Newtonian hydrodynamics, time and space are understood as separate entities. Only one 
definition of the time slices is possible: orthogonal to the time axis. However, a free choice of 
coordinates within the time slices persists. We may choose Eulerian coordinates along the time axis 
or Lagrangian coordinates along fiuid element trajectories. If one understands the grid point labels 
in a dynamical adaptive grid technique as coordinates, many mappings between such coordinates 
and the location of fluid elements can be realized with appropriate prescriptions for the motion of 
the adaptive grid. These choices are easily interpreted in the general relativistic view: The Eulerian 
coordinates correspond to vanishing shift vectors and the Lagrangian coordinates correspond to shift 
vectors equal to the three-velocity of fluid elements. In the adaptive grid case, the shift vectors are 
identified with the grid velocity. 
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In the following investigation, the time slices are chosen orthogonal to matter trajectories such 
as to produce comoving coordinates with vanishing shift vectors. Nevertheless, the relationship 
between adaptive grids in Newtonian space-time and shift vectors in the general relativistic context 
provides a guideline for the implementation of shift vectors in general relativistic simulations. Fixed 
coordinate intervals can be mapped to the time slices at will in order to resolve physically interesting 
regions with high precision. 

We start with the spherically symmetric Einstein equations in comoving orthogonal coordinates 
as given by Misner & Sharp (1964). They are based on the metric 



ds'^ = -a^dt^ + (^-^ da^ + {d^'^ + sin^ M^p"^) , 



(2) 



where r is the areal radius and a is a label corresponding to an enclosed rest mass (the prime 
denotes a derivative with respect to a: r' = Or /da). The proper time lapse of an observer attached 
to the motion of rest mass is related to the coordinate time dt by the lapse function a. The angles 
■& and If describe a two-sphere. 

The thermodynamic state of a fluid element is given in its rest frame by the rest mass density, 
p, the temperature, T, and, in our application, the electron fraction Y^. An equation of state renders 
the composition in nuclear statistical equilibrium, the specific internal energy, e, and the isotropic 
pressure p. In analogy to the definitions of Romero et al. (1996) in polar slicing and radial gauge, 
we compose conserved quantities for specific volume, total energy, and radial momentum: 

r = Te+ - — - l-u^ (4) 



r + 1 

S = u{l + e). (5) 

The velocity u is equivalent to the r -component of the fluid four- velocity as observed from a frame 
at constant areal radius (May Sz White 1967). In the special relativistic limit, T = ^yl + u'^ — 2m/r 
becomes the Lorentz factor that takes the boost between the inertial and the comoving observers 
into account. The gravitational mass m is the total energy enclosed in the sphere at rest mass a. 
In the nonrelativistic limit, we retrieve with a = T = 1 the familiar specific volume 1/D = 1/p, the 
sum of the specific internal, kinetic, and gravitational energy t = e + — m/r, and the specific 
radial momentum S = u. 

With these definitions, the complete system of general relativistic hydrodynamics equations 
can be written in a conservative form (Liebendorfer, Mezzacappa, & Thielemann 2001): 



d_ 
di 



1 

D 



d 

— — kvrr^aul (6) 
oa ^ 



- = --[A^raup] (7) 



-8- 



d 



da 



a 
r 
1 

D 



3p 

l + e + — 
P 



— + 8Trr^(l + e)p- — 
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' = ^'^'^Ta^-pd-a^^'^- 



(8) 
(9) 
(10) 

(11) 



Equations (6)- (8) describe the conservation of volume (analogously to mass conservation in Eulerian 
coordinates), total energy, and radial momentum respectively. The constraints (9) and (10) explain 
themselves in the analogy to the Newtonian limit, where the first becomes the definition of the rest 
mass density and the second the Poisson equation for the gravitational potential. The enclosed 
volume in Eq. (9) is defined by the areal radius as in the Newtonian limit, V = Airr^/S. The 
condition for the lapse function in equation (11) is derived from the space-space components of the 
Einstein field equations. 



2.2. Shift vectors and the dynamical adaptive grid 

In Fig. (1) we have now moved from the equations of general relativistic hydrodynamics to 
their conservative formulation in specific coordinates as shown in the lower right corner. In order to 
traverse to the left towards a stable implementation of the adaptive grid, we rewrite the derivations 
of Winkler, Norman, &; Mihalas (1984) in a form that applies to the general relativistic case as 
well. The goal is a stable and conservative finite difference representation of equations (6)- (8) in 
terms of n grid points, {ai{t),i = 1 . . . n), that continuously move with respect to the enclosed rest 
mass label. 

We concentrate on observers at rest in their slice. They reside in local orthogonal reference 
frames that we assume to be coUinear with global coordinates A. Consequently, coordinates A 
have vanishing shift vectors everywhere. We introduce another system of global coordinates, B, 
having arbitrary shift vectors on the same time slices. Let us select a fixed grid in coordinates B: 
{g^gN3 G R-^) H = 1 . . . ni, Z2 = 1 ... 712, ^3 = 1 • • • ns}. The trace of the world-lines along this grid 
defines a moving grid in coordinate system A: {a{t,qi)}. The grid world lines along coordinate 
system B cut the continuum into zones Aq on each time slice. These zones contain a time dependent 
selection of observers in coordinate system A: Aa = {a{t,q),q G Aq}. Based on the observations 
y of single observers a in the zone, we define a zone integral of the observations 

{yit,d))= [ y{t,a)da. (12) 

JAa 

In coordinate system A, the borders of the zones change with time according to the shift vectors 
in coordinate system B. In Newtonian parlance, one can identify the moving zones in coordinate 
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system A with the motion of an adaptive grid. The important point is that the adaptive grid only 
regroups the observers into new zones. As stated akeady by Winkler, Norman, & Mihalas (1984), 
the adaptive grid never changes the reference frame of the observations. It solely determines the 
zonewise support for the integration of the observations in the time slice. 

In order to formulate time evolution equations, we are interested in the temporal change of 
the zone-integrated quantity 

{yit,q)) = [ y{t,a)da (13) 

in a specific zone lS.q. We can apply the Reynolds theorem to relate the time derivative d{y(t, q)) /dt 
in coordinate system B to the time derivative {dy{t,a)/dt) in coordinate system A: 



d{yit,q)) 
dt 



/ dy{t,, 
\ dt 



+ 



y{t,d) 



d{Aa) 



dd{t, q) 
dt 



n 



dS. 



(14) 



Note that the Reynolds-theorem is a mathematical relation that has no physical content. The left 
hand term is the temporal change in the zone average of y in the numerically accessible coordinate 
system B. The first term on the right hand side relates this time derivative to the time derivative 
in orthogonal coordinates A that is easily evaluated based on the physical equations describing 
the time evolution. The correction in the second term on the right hand side is due to the grid 
motion. It is an integral over the zone surface elements dS (with normal n) that corrects for the 
observations made by observers that enter or leave the zone during grid motion. 



If an equation describing the physical evolution is written in a conservative form 



/ dy{t, a) 
\ dt 



+ 



9(Aa) 



/(t, a) ■ fi \ dS = I s{t, d)dd 

^ JAa 



with a source s{t,a), its formulation on the adaptive grid simply becomes 
d{y{t,q)) 



dt 



+ 



9(Aa) 



j^, . ^.dd{t,q) 

fit, a) - y{t, a) 



■ n 



dS 



/ s{t, a)da. 

JAa 



(15) 



(16) 



The motion of zone boundaries naturally leads to a fiux entering or leaving the zone. The rate 
of coordinate change at a fixed grid point v^^^ = — defines a "grid velocity" with respect to 

coordinates A. In our spherically symmetric case, the Reynolds theorem reduces to 



d_ 
dt 



a-{t,qi+i) 



j'a(t,qi- 

Ja{t,qi) 



yda 



a(t,qi) 



dt 



da + 



d a{t,q) 

dt 



y- 



a{t,qi. 



- a(t,qi) 



and leads for example to the continuity equation (6) on the adaptive grid 



g ra(t,qi+i) J 
9i Ja(t,q,) D 



da + 



2 1 da{t,q) 

—4:Trr au — — — 

D dt 



a(t,qi) 



= 0. 



(17) 



(18) 



The integral of the conserved quantity (e.g. volume) over the rest mass of the star only depends 
on surface terms in this favorable form because contributions from zone interfaces exhibit exact 
cancellation. Thus, the motion of the adaptive grid affects local resolution but not global conser- 
vation. 
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2.3. Dynamics of the adaptive grid 

After having described the formulation of physical equations on an adaptive grid, we focus on 
the independent question, how to adapt the grid to interesting features in the physical evolution. 
The basic idea is not only to resolve steep gradients with a locally increased grid point concentration. 
To a certain extent, the dynamical adaptive grid is also capable to captTirc a self-similar flow and 
propagate it through the computational domain mainly by grid motion instead of rapid changes in 
the physical state of zones. This can lead to larger time steps when steep gradients are present. 
For example, if a shock front is smeared out over several grid points, it takes a width Aq in grid 
point labels q (see Fig. (2)). Suppose the shock moves with velocity Vq with respect to the grid; 
then there is time At = Aq/vq to change the variables from their pre-shock to their post-shock 
values. We consider the variable y with the largest relative jump \dy/y\ at the shock front. If the 
relative change of the variable per time step is limited to e one gets from + < \y\[l + e)"^ 
about n = ln(l + \dy / y\) / Xivil + e) time steps for this change. Therefore, the time step is limited 

^* < f'''°"ir',v (19) 

In contrast to Lagrangian or Eulerian schemes, this time limit is not restrictive with the adaptive 
grid. Firstly, the high resolution at the shock front gives a large Aq and secondly, the fact that the 
grid points tend to move with the shock structure makes the velocity difference between the shock 
front and the grid, Vq, very small. 

We have implemented the grid equation of Dorfi &: Drury (1987). We only differ in one detail: 
The scheme has been transferred from Eulerian coordinates in flat space to comoving coordinates 
in Misner- Sharp time slices. Thus, we retrieve a Lagrangian scheme if the adaptive grid is switched 
off. The grid equation of Dorfi and Drury is based on a resolution function, 




where the sum includes a selection of relevant variables y as e.g. velocity, density, etc.. We set 
the global scale in this equation to y^'^^ = \yi+i ~ yi\- In terms of a grid point concentration 
n = a^'^^{da{t, q) /dq)~^ the grid equation implements the requirement 

If we imagine the variables y/y^'^^ plotted in a graph versus a/a^'^^, this condition enforces a constant 
value of 



— - 

n _ ( f da{t,q)Y f dy{t,qy 



R \\ aP^^dq J ^ \ y^^dq 

equivalent to a constant generalized arc length per grid point interval Aq. The projection of grid 
points to the a/a^'^'-axis results in an increased grid point concentration at locations with steep 
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Fig. 2. — This shock tube simulation (described in section 4.5) illustrates the dynamics of the 
adaptive grid. Data at two different times are shown. Graph (a) shows the scaled velocity profile 
in coordinates A, i.e. versus enclosed mass. We indicate the velocity jump across the shock with 
dy. The grid is set to produce equidistant generalized arclength between grid points. This is only 
qualitatively visible in this graph because of the influence of the other variables and the spatial 
smoothing operator that guarantees a smooth reduction of the resolution with increasing distance 
from the shock. Graph (b) shows the same profile in coordinates B, i.e. versus grid label. The 
shock is smoothly spread over an interval Aq that is used to estimate the time step restriction, 
equation (19). The profile in the computer representation (b) is fairly stationary, while the shock 
propagation in physical space (a) is dominated by the evolution of the mapping a{t,q). 
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gradients in the variables y/y . This basic approach is refined by the sophisticated apphcation 
of operators for spatial smoothing and temporal retardation of the grid adaption as described by 
Dorfi & Drury (1987) (see Fig. (2)). 

3. Numerical Method 

The nonlinear partial differential equations (6)-(8) and the constraints (9)-(ll), (21) are finite 
differenced in the spatial direction to build a large system of coupled ordinary differential equations 
with algebraic constraints. Each independent variable in each zone acts as an individual unkown 
function of time in the solution vector. In a second step, this large system of potentially stiff 
algebro-differential equations is integrated in time with our solver AGILE. Experimentation was 
necessary to find the final finite difference representation of the hydro dynamical equations. We 
have therefore separated the solver for the ordinary differential equations and the adaptive grid 
equation from the modules that describe and discretize the actual physical equations. AGILE has 
the capapility to calculate the Jacobian of the physical equations automatically for the implicit 
time evolution. This feature substitutes the cumbersome and errorprone coding of the Jacobian 
and leaves room for the motivation to change and improve the physical equations. 

In the following subsection we start with the description of the time evolution in AGILE. 
Then, we discuss the numerical diffusion induced by the advection that accompanies the adaptive 
grid. We introduce artificial viscosity and review the role it plays in combination with the adaptive 
grid technique. Finally, we outline our complete spatial discretization of the Einstein equations in 
comoving coordinates. 

3.1. Time evolution in AGILE 

We start with a vector, y, that depends on the parameter time, t. The temporal evolution of 
the vector y shall be described by an implicit system of equations 

f(<,.,|)=0. (23) 

The set of equations F is assumed to contain the same number of equations as the vector y contains 
components to be solved for. In the discretization of these equations we use upper indices to denote 
the time steps: is a vector of the variables at time and j/""*"^ is the vector after the time step 
dt = t'^'^^ — f^. The discretized equations then have the form F (y", y"''"-'^; dt) = 0, where the time 
derivative of the state vector y is calculated by a finite difference representation included in the 
equations F. 

A time step dt is chosen such that the relative change in the components of the state vector 
is not expected to exceed a given value in a single time step. The solution vector at time 
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f-n+i _ _|_ jg ii^Qj^ found by a Newton-Raphson scheme: The equations are Taylor expanded 
around a guessed solution vector 

+ Ay; dt) = F{y\ dt) + ^^^^^^^^^Ay + 0{{Ayf). (24) 

and solved for the corrections 

A»^-(^^^^9P^)"^fe».^-;*). (^5) 

This procedure is iterated until a norm of the correction vector \/\y\/y^'^'^ becomes sufficiently 
small. We use the scaling y"'^"^! = |?/"| + The first term usually dominates and the second 

term limits the sensitivity to numerical noise in a very small component. 

The chosen notation does not show at which time between and t""*"^ the equations actually 
are solved. In the fully implicit approach, the old values are only used in the time derivatives 
and the equations are solved at f""*"^. In the literature one also finds semi-implicit approaches 
(Winkler & Norman 1986; Mair 1990; Swcsty 1995) that solve the equations at a time t"^'^ 
depending on some averaged variables y"+-^(y"^ y"+i) with < A < 1. For A = one obviously 
gets an explicit scheme. Optimally, with A = 1/2, a semi-implicit scheme is achieved whose time 
discretization error is second order in the time step. Unfortunately this choice is often not stable 
enough. A compromise are schemes with A ~ 0.7 which is reported to be about the minimum A with 
satisfying stability (Winkler &; Norman 1986). For the time being, we use the fully implicit choice 
A = 1 because this setting is sTifficient for the supernova application. However, we experimented 
with a semi-implicit multi-step extrapolation method that generalized the approach of Bader & 
Deufihard (1983) to algebro-differential equations (Liebendorfer & Rosswog 1998; Liebendorfer 
2000). The great benefit from the higher order extrapolation in combination with the smooth 
computer representation of the physical state on the adaptive grid was, however, severely disturbed 
by the discontinuous switches in the upwind advcction scheme. We therefore focus on the simple 
and stable fully implicit integration and revive our higher order integration scheme in AGILE 
(Adaptive Grid with Imlicit Leap Extrapolation) at a later occasion. 

We implemented the following automatic procedure that calculates the Jacobian numerically 
based on the modules containing the physical equations. For the first iteration, we assume that 
we have a guess y""*"^ (e.g. y""*"^ = y'^) at hand for the state vector at the new time and that 
we know the sparsity pattern of the Jacobian well enough to specify a region that contains all 
nonzero coefficients (e.g. a large band). At the old time we then attribute a scale y"'^'^' to the 
state vector y" and evaluate the residuum vector R = F{y"-, y"'^^;dt) based on the guess. For the 
calculation of the numerical derivatives, we define a small number e based on machine precision. 
Furthermore, the components of the state vector, denoted with label i, are sorted into J distinct 
groups gj according to the rule that none of the equations F may depend on more than one state 
vector component out of the same group. Based on the sparsity pattern, the groups are formed 
such as to minimize the required number of distinct groups. We select a group gj and create a copy 
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yvar,n+i ^ ^n+i q£ guess with a Small variation 



A[z] = (26) 

in all components i G gj. The residuum vector K""'^ = F{y"',y""'^'"'^^;dt) is then calculated with 
the varied guess. Prom the two residua we can extract the components of the Jacobian 

^[^A-y[A - ^[.j (27) 

for all i ^ Qj. The Jacobian is complete when this procedure has been performed for all groups J. 
Finally, we scale the rows of A and the right hand side residuum vector by the maximum component 
in the row and solve the linear system 

A[K i] (Ay[i]/y[i]"'-') = -R[k]. (28) 

i 

to get the corrections to the guess. Before the next iteration, all zeros in the Jacobian are detected 
and the sparsity structure is refined. This allows further reduction of the number of required groups 
for all following time steps. The evaluation of By"-'^ = F{y"', dt) for a group j is completely 

independent from the corresponding calculation for a different group j' ^ j. In order to compose 
a complete Jacobian, the system of equations has to be calculated once with the actual guess, and 
J times in parallel with varied guesses. 

This parallelism, in principle, allows the build of the Jacobian on parallel machines without 

affecting the modules with the physical equations. The whole system of equations is still evaluated 
at once. The physical equations for the different components in the state vector can therefore make 
use of common subexpressions or be vectorized. If for each group a separate process is available, 
the described parallelism maximally reduces the wall-clock time per implicit time step to the single 
process time of the evaluation for one residuum vector as in an explicit scheme. However, this 
comparison does not include the time spent with the inversion of the Jacobian. A realistic estimate 
of the gain in wall clock time has to await a physical need to complete the implementation of 
this feature. Note also that the described separation into groups is not efficient with an arbitrary 
sparse structure. It however is efficient with the block-diagonal Jacobian occurring in spherically 
symmetric hydrodynamics. 



3.2. Diffusion and artificial viscosity 

Before describing a detailed finite difference representation of the Einstein equations (6)-(ll) 
in the next section, we stress here a few discretization principles that were useful in supernova 
simulations, and might carry over to other applications as well. The modeling of core collapse 
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supernovae presents a special challenge for the finite differencing of the hydrodynamics equations. 
On the one hand, global energy conservation is important because of the comparable size of the 
gravitational, internal, and neutrino radiation energy. The balance between those defines the two 
orders of magnitude smaller total energy. The latter is comparable to the expected explosion energy 
and has to be accurately evolved in order to allow conclusions on the explosion energy in the model. 
On the other hand, the accurate evolution of the internal energy and temperature is essential for 
the strongly temperature-sensitive weak interaction cross sections that figure as a prerequisit to 
neutrino transport, which may determine the outcome of a supernova simulation. Thus, one would 
like to solve an equation for the evolution of global energy to obtain exact energy conservation, and 
an equation for the evolution of internal energy to obtain accurate temperatures, and a momentum 
equation to determine accurate fiuid velocities that relate to the kinetic energy. However, this set of 
equations is overdetermined since the total energy (without gravitation) is nothing more than the 
sum of internal and kinetic energy. The following measures were necessary to obtain good results 
in all three quantities with the solution of only two independent equations: 

(i) In the Newtonian limit, the time derivative of the specific total energy, r = e + l/2u'^, 
relates its evolution to the evolution of specific internal energy, e, and the evolution of specific 
linear momentum (=velocity), u. Consequently, the finite difference representation of the internal 
energy equation and the momentum equation determine a consistent finite difference representation 
for the total energy equation. The challenge is to find a discretization of these constituents such 
that energy conservation becomes manifest in the finite difference equation for total energy. By 
making the three evolution equations consistent in their finite difference representation, it is less 
important, which two equations are independently evolved and which third quantity is determined 
by them. 

(ii) However, we found an exact match only for the terms up to order (v/c). Moreover, the 
match does not include the effect of the advection terms, and truncation errors may swamp the 
solution of the dependent quantity when it happens to be generated out of the cancellation of two 
much larger independently evolved quantities (e.g. the evolution of total energy and large kinetic 
energy lead to a poor extraction of internal energy). In order to keep these inaccuracies small, we 
pick the evolution of total energy as our first independent equation and complement it with the 
evolution of the ratio of kinetic and internal energy as second independent equation. The time 
derivative of this ratio, 



evolution equation of the smaller quantity. We consider global energy conservation and accurate 
internal energies to be more important in a supernova simulation than the lost strict conservation 
of radial momentum. 




(29) 



is a linear combination of the internal energy equation and the momentum equation that favors the 



(iii) The finite difference representation of the dynamical equations introduces numerical dif- 
fusion. Part of this diffusion is desired to stabilize the solution. E.g. the upwind differencing in the 
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momentum advection scheme acts as artificial diffusion, dissipates kinetic energy, and consequently 
leads to a mismatch in the relation between total, internal, and kinetic energy evolution. However, 
the artificial diffusion introduced by upwind differencing can be evaluated and related to a source 
term for the internal energy equation. The dissipated kinetic energy is then correctly transformed 
to internal energy that restores overall consistency. We will include this diffusive advection term 
into the concept of artificial viscosity that is needed in our scheme to limit the resolution of a shock 
front. The artificial viscosity tensor is based on physical viscosity and consistently included in the 
general relativistic set of hydrodynamics equations. It provides the desired mechanism for dissipa- 
tion of kinetic energy into internal energy. A more detailed discussion of the artificial diffusion and 
viscosity is the aim of this section. 

Along the lines of equation (17), a generic time evolution equation for a conserved quantity y 
would be finite differenced as 



Vi'dai' - Vi'daj' ^dv „,adv „,ext 

It +^^+1 



Here we denote the spatial dependence of the quantities by zone edge indices i and adress a zone 
center value by the index i' = i + 1/2. Quantities with an overbar belong to the past time t while 
all other quantities belong to the present time t + dt. The temporal change of the quantity yi'dai' 
in a zone with extent doj/ is given by the advection at its boundaries and a not further specified 
external source y^^^. Let us set the advection terms by first order upwind differencing depending 
on the direction of the relative velocity u^^^ defined in equation (41): 

* I yi'uY^ otherwise. 



We would obtain an identical advection scheme by writing 

(yi' - yi'-i) . (32) 



yf^ = \uf{y,+y,^,)-\ 



da. (33) 



Hence, we may interprete it as the superposition of the second order accurate interpolated advection 
'^T^ {yv + yv-i) /2 plus a correction term — \uf^\ {yi' — yi'-i) /2. It is well-known that the second 
order term alone does not lead to a monotonic and stable advection scheme. The additional 
stabilizing correction term is a diffusive flux that corresponds to a diffusivity 

1 
2 

We therewith note that the advection (related to the dynamical adaptive grid) introduces a diffu- 
sivity that is proportional to the velocity of the grid with respect to matter, and proportional to 
the grid spacing. In the case of the momentum equation, this diffusivity reduces the kinetic energy 
and consequently leads to internal energy generation if conservation of total energy is supposed to 
hold. In the adiabatic Newtonian limit, only the internal energy equation 



de ( 1 



5t+ ^+2 



A.||(i|=o (34) 
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complements the momentum equation consistently. In the finite difference expression of the general 
relativistic momentum equation we take this effect into account by not directly coding upwind 
differences. Almost equivalently, we use the second order advection term in equation (46) and 
add the stabilizing diffusion term (33) to the artificial viscosity tensor (48). Its role as source for 
internal energy is then automatically accounted for by the entry of viscous pressure in the energy 
equations. The implementation conceptually mutates from a diffusive implementation of ideal fiuid 
dynamics to an "ideal" implementation of diffusive fiuid dynamics. 

Artificial viscosity is needed to handle discontinuities in the solution that may occur in ideal 
fiuid dynamics, e.g. in shock waves. In a simple finite difference scheme on an equidistant grid, 
numerical oscillations appear as soon as the steepness of the shock front exceeds the resolution of 
the grid (Noh 1978). These oscillations can be suppressed by the introduction of artificial viscosity 
that limits the steepness of the shock wave to the resolution that is provided by the grid (Von 
Neumann k, Richtmyer 1950). The adaptive grid puts the same concept into a slightly different 
light: The grid adapts according to the steepening gradient and is always able to resolve it properly 
without initially leading to unphysical oscillations. However, with increasing resolution, the mass in 
the zones at the shock front approaches zero and the system of equations becomes ill-conditioned. 
Poor convergence of the implicit step is the immediate response. Artificial viscosity is needed in 
our context to limit the gradient at the shock front and prevent the grid resolution from growing 
indefinitely. 

In spherical symmetry, however, the formulation of artificial viscosity has to be chosen with 
care to circumvent systematic artificial heating in phases of homologous compression. We extend 
the approach of Tscharnuter k, Winkler (1979) to general relativistic hydrodynamics (Liebendorfer, 
Mezzacappa, & Thielemann 2001) and define the viscous pressure based on a parameter AZ with 
the dimension of a shock width, 

Q ^ [ A/Vdiv(«) [f - idiv(«)] ff div(u) < 0, 

I otherwise 

d 

div(it) = (47rr^«) . 

The corresponding viscous heating in the energy equation becomes: 

de _ f u du\ Q 
adt \r dr J p 



du 1 . ^ 



^. (36) 
P 



The first expression shows that the viscous heating vanishes in the case of homologous compression 
(u/r = du/dr). The second expression together with equation (35) shows that viscous heating 
always has positive sign. In our system of equations, the viscosity affects the equation for the total 
energy evolution (7), the momentum evolution (8), and the constraint for the lapse function (11). 
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These equations become (Liebendorfer, Mezzacappa, & Thielemann 2001) 

[A-rrr^au {p + Q)] 
[Anr'^ar {p + Q)] 



dr 
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at 
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!!l + 87rr2(l + e)(;, + g)-^ + ^ 
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J_d_ 
Vp aa 



[VaQ] . 



(37) 



3.3. Discretization of the Einstein equations 

In all the other sections, the velocity of light, c, and the gravitational constant, G, were 
absorbed in the choice of the units. For clarity, we include these constants explicitly in this section, 
where the details of the finite differencing are given. As before, we denote by i the zone edges 
and by i' = i + 1/2 the zone centers. Quantities that are evaluated at the old time are marked 
with an overbar, all other quantities belong to the new time f^^^ = + dt. The state vector y = 
{a,r,u,m, p,T,Ye,a} contains the independent variables: enclosed baryon mass, radius, velocity, 
enclosed gravitational mass, rest mass density, temperature, electron fraction, and lapse function, 
respectively. Prom the equation of state we get the pressure p and specific internal energy e. First, 
we specify the enclosed volume and its time derivative 

Vi = -j^, Wi = Anrfui, (38) 
and define zone differences for the radius and enclosed rest mass 



dvi' = n+i - Vi, daii = Oj+i - Oj. (39) 

Further, we define = ^/l + [uijcf - 2Gmi/{c^ri) on zone edges and prepare the conserved 
"density," and the "radial momentum," Si. Wc split the conserved energy into "internal" 
energy, rlj/, at the zone centers and "kinetic" and "gravitational" energy, r2j, T3i, at the zone 
edges 



Di, = 



Pi' 
Ti' 

= (^ + 5)"^ 



rl,;' = Fj/e 
r2,: 



2 



Vi + l 

Fi + 1 ^ ' 
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We interpolate F to the zone centers and e to zone edges by arithmetic means. Based on the deriva- 
tion in section 2.2 and equation (17), the adaptive grid causes advection of conserved quantities 
through the zone boundaries. We denote the relative velocity between zone edges and matter with 

= (41) 

and implement first order upwind differencing for the advection terms. On the zone edges this leads 
to 



lef = Y^^^ruf (42) 



if uf' > and to 



A' ^ 

Tlf- = Tl,uf 

Yif- = Ye,i>v^^' (43) 
otherwise. At the zone centers we set u^f = (n^^^ + "^T^)/^ ^^'^ advect the quantities 

T2f'' = T2iuf (44) 

if uf > and 

T2f^ = T2i+^uf (45) 
otherwise. For the advection of gravitational energy and momentum, we set 

Sf^ = ^iS,+i + S,)v^\ (46) 

Since the gravitational energy is a rather smooth function, that depends only weakly on the im- 
mediate solution of the energy and momentum equation, we can implement the more accurate, 
direction independent, second order advection term. The radial momentum advection, however, 
is not smooth. But we prefer to split the upwind differencing into second order advection and 
a diffusive term that is included in the artificial viscosity. The latter is based on the velocity 
divergence 

4:7T 

dwu, = min (o, (47) 
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and is given by equation (35) and (33) 



Qi> = Afpi/divui' ( ^^^^^.^ ~ -^diviij/ ) (48) 



(Si+i — 



aiiTiidVi' 



Q 3 / Ui+i -Ui 1 \ Qj/ 

^i' = ~o — ^ odivMi' dai 

2 \ dri> 3 y ft' 



After these preparations, we begin with the finite diff'erencing of the constraints from equations 
(9)-(ll). Straightforward is the discretization of the definition of the density and the Poisson 
equation 



dtti' - Di'dVi' = (49) 

2 



nii+i — + Tj/ (1 + ^ ) dai' = 0. (50) 



Less straightforward is the implementation of the constraint for the lapse function that is in some 
details arbitrarily chosen to become 

ai'Pi' - ai'^ipi'-i + — {Vi'ai'Qi' - Vi'-iai>-iQi'-i) 

-^a^'.i^r* + ftc^ (l + ^) (a,, - a,,_i) = 0, (51) 

where the zone center values of V and the zone edge values of the density p and internal energy e 
are found by arithmetic means. 

We can then start with the construction of the time evolution equations. The continuity 
equation lives on zone centers and is discretized as 

+ Vt^l - Vr"-" - {ai+iwi+i - aiWi) = 0. (52) 

The equation for internal energy evolution also lives on zone centers and reads 

_ ei'dttj' -Ci'daj' ^ 

+ ^{wi+^-Wi)-ai,e'^,-ai,ef' = Q. (53) 

t i' 

We introduced an additional energy source e?,^* that describes energy exchange with external pro- 
cesses as e.g. with a nuclear reaction network or a radiation field. Analogously, we also include an 
external stress Sf^^ and an external compositional change Y^^} . The momentum equation lives on 
zone edges and is discretized as 

jpS _ Sjdai — Sjdai 

— r ^i' ~ '^i'-i 



- 21 - 



+ — [Vi {Ti'ai'Pi' - rj'_ia 



+ 
+ 

+ 



{TiiVi>ai'Q. 
ai \ / e 
n \ c 

.2 



ipi'-l) 

^i'-iVi'-iaii-iQit-i)] 
QVi {Pi + Qi)\ Grrii 



1 + 



1 



ui {2pi - Qi) - {pi + Qi 



(54) 



We now use the sum of the internal energy equation and velocity times the momentum equation 
in the Newtonian limit as a guideline for the discretization of the total energy equation. In the 
Newtonian limit, and with the omission of external sources, advection, and artificial viscosity we 
find 



Ei'dai' — Ci'dai' 
Jt 

Ui+idui+i - 



+ Pi' {wi+i - Wi) 
- Ui+idai+i 



dt 



Ui+iGnti+i 
+ m+i [Pi'+i -Pi') + dai+i. 



(55) 



The expression Wi^iPi' from the internal energy equation cancels with its counterpart from the 
momentum equation and leaves the adiabatic work Wi+iPi'^i — WiPi' that is applied to the fluid 
element. Which fluid element? We observe that the conservation of energy is expressed by relating 
the internal energy enclosed by zone edges, and the kinetic and gravitational energy enclosed by 
zone centers, to this adiabatic work. The work excerted on such a "distorted zone" is given by the 
product of the vohmic changes w at zone edges with the pressure p at zone centers. The position 
of these terms perfectly fits the arrangement of the variables on the staggered grid (and has to 
be remembered when one implements energy conservation checks). This suggests the following 
discretization of the full energy equation (7) 



{rliidaii + T2i+idai+i + rSi+idui+i) — 

- (rlj/daj' + r2j+i(iai+i + r3j+idaj+i) ^ 



-Wiai' 



2 (Pi' +Pi') + Qi 



- {pi>+i +Pi'+i) +Qi'+i 



„, -p „ext 



0. 



(56) 



It is problematic to apply time centering on the adaptive grid because the quantities belong to 
different fluid elements across a time step. The finite differencing is therefore kept first order 
backward Euler. Two exceptions, however, improved the accuracy of the solution considerably 
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while they were not found to conflict with the grid motion. In the total energy equation, we time- 
average the pressure as in (Swesty 1995) and in the momentum equation we time-average the 
term in the gravitational force as in (Yamada 1997). We solve the total energy equation (56) 
together with the momentum equation (54) or the internal energy equation (53) depending on the 
contribution of internal and kinetic energy to the total energy. As discussed in the previous section, 
the momentum equation and the internal energy equation are merged with weights according to 

e^'F^il - \Si+iFl, = 0. (57) 

The evolution of the electron fraction is determined by that describes weak interactions with 

neutrinos 

Yp jidcLji — Y p jidoLi' 



' ' + Y,il^ - - cc,Y!f = 0. (58) 

The eight equations (49)-(52) and (56)-(58) are implicitly solved in one time step, together with 
the adaptive grid equation (21). 



3.4. Boundary conditions 

For completeness, we add in this section details about the implementation of the boundary 
conditions at the border of the computational domain. The boundary conditions are partially 
modified for the successful run of the individual test problems as described below. 

We set the zone edge with i = 1 at the center of spherical symmetry and impose boundary 
conditions 

ai-ai = (59) 
n - f 1 = (60) 
ui — ui = (61) 
m\ — nil = (62) 

for the variables that live on zone edges. The reason for rather keeping them constant instead of 
explicitly setting them to zero is that this formulation remains valid in simulations where the inner 
boundary is off center, as e.g. in the shock tube calculations presented in the next section. Special 
care has to be taken of the advective fluxes at the zone center z = 1 -|- 1/2. Although the kinetic 
and potential energy, as well as the radial momentum, vanish at the center, the content of the first 
half zone with respect to these conserved quantities is not zero. It is subject to advection if the 
zone center i = 1 + 1/2 moves inwards. In the following reasoning about the distribution of these 
quantities around the symmetry center, we will neglect the difference between the gravitational 
and rest mass, and set T = 1. We assume that a/r^ = A is constant in the neighborhood of the 
symmetry center and estimate for the gravitational energy in a sphere 

-da = -A3a3. (63) 
r 5 



L 
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For the evaluation of the kinetic energy and radial momentum, we assume u/r = B constant and 
evaluate 

1 9 , 3 „9 ^ 2 5 

-u da = -T7:B A sas 
2 10 

3 14 

uda = -BA~3a3. (64) 
4 



Jo 



The amount of the conserved quantity in the first half zone is then found by integrating the 
conserved quantity over the first one and a half zones from ai = to 02' and subtracting the 
content that has already been taken account of in the range da2- With an equidistant zoning, 
= 3/2da2, we find the relations 



'"'2' a , 3 /3\ t 02 



Jo 



/ —da = - ( - ) —dao 
Jo r 5 \2j r2 

1 2, 3 /3\t 1 2 , 

-u da = - - 7.U2da2 



f 

Jo 



2 5 \2J 2 

"2' 3 /3\3 

uda = ^ ( 2 ) '^2da2 (65) 



that express the content of the innermost one and a half zones in terms of the content in the interval 
da2 that is well placed on the staggered grid. We therefore define the fractions /e = 3/5(3/2)^/^ — 1, 
fs = 3/4(3/2)^/^ — 1, and set for advection purposes a content 

t2i = /eT22 
t3i = /eT32 

Si = fsS2 (66) 

for the kinetic energy, gravitational energy, and the radial momentum in the innermost half zone. 
These measures are not required in the simulations of the shock tubes, where the inner boundary 
is off center and the physical state constant in its neighborhood. 

At the outer boundary, we impose the conditions 

(67) 
(68) 
(69) 

«n'-(i-^)r;/ = (70) 

(71) 

We determine the quantities m, r, and u on the last zone center n' from the last zone edge n with 
the equations (49), (50), and Wn' — Wn = 0- Equations (67)-(69) implement a constant surface 



Pn' 


- Pn' 


= 


Tn' - 


- Tn' 


= 




^ e,n' 


= 


2m„' ^ 
rn' / 


) 


= 




-an 


= 



-24- 



pressure boundary condition by the equation of state. Equation (70) matches the lapse function to 
the exterior Schwarzschild metric. Note that in equations (59)-(62) and (67)-(71) every independent 
variable fulfills one boundary condition, except for the rest mass a. The computational domain 
is tied to the desired mass range by two boundary conditions that compensate for the missing 
boundary condition for n/R in the adaptive grid equation (21). The constant pressure boundary 
condition is not very practical for the homologous collapse test because the surface tends to drift 
outwards during the time evolution and destroys the homologous velocity profile. We therefore 
evolve the surface pressure proportional to the change in the central pressure and replace equation 
(67) and (68) by 

Pn'-^P7' = (72) 
Pi' 

Tn'/pl7'-T^'/py = 0. (73) 
The latter condition imposes constant entropy at the surface. 

4. Supernova-Related Standard Test Problems 

In this section we demonstrate the capability of our hydrodynamics code to accurately repro- 
duce important features in stellar core collapse and supernova explosions. We present the standard 
test calculations for supernova codes (Swesty 1995; Yamada 1997) that complement the appli- 
cation of AGILE in (Mezzacappa et al. 2001) and (Liebendorfer et al. 2001). All examples are 
based on a polytropic equation of state and, unless marked otherwise, simulated with the same 
general relativistic code version. 

4.1. Sedov point-blast explosion 

A nonrelativistic point-blast explosion in spherical symmetry was analytically analyzed by 
Taylor (1950) and Sedov (1959) (we followed Landau & Lifschitz (1991) for the reproduction 
of the self-similar analytical solution to this problem). An amount of energy Etot is deposited in 
a uniform gas with a density po and internal energy that is negligible with respect to Etot- We 
start with an initial specific energy cq = 10~^ at a density po = 1 in cgs units. The pressure is 
given by the ideal gas equation of state p = (7 — l)pe with 7 = 5/3. The computational domain 
has a radius of r = 1 around the center where the energy Etot = 1 is deposited. The initial state 
is prepared by starting with 100 equidistant grid points. The explosion energy Etot is placed into 
the innermost zones according to an exponentially decaying spatial distribution on a time scale 
between the dynamical time scale and the much shorter time scale of the grid adaption. By this 
measure, the adaptive grid can move inwards and resolve the exponential shape of the deposited 
energy. If the amount of Etot is resolved, the time scale increases to the dynamical time scale 
and the shock wave starts to move radially outwards. The numerical solution is compared to the 
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Fig. 3. — This figure shows the Sedov point blast explosion, the spherically symmetric explosion 
after the deposition of a large energy amount in the symmetry center at r = 0. The dots indicate 
the location of the grid points in the numerical solution. They are compared to the solid line, which 
represents the analytical solution. 
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analytical solution at time t = 0.5 in Fig. (3). We find the correct shock speed and emphasize the 
resolution achieved with the adaptive grid. The deviations in the inner part of the sphere stem 
from the nonideal initial conditions. 



4.2. Nonrelativistic and relativistic shock tube 



Probabely the most popular test is the calculation of a shock tube. The initial conditions for 
a nonrelativistic simulation were proposed by Sod (1978) and a special relativistic version was 
calculated by Centrella & Wilson (1984) and Marti & Miiller (1994). Sod's shock tube evolves 
the decay of a discontinuity in density p and specific energy e in a closed box. The box is filled 
with an ideal gas whose pressure is given by a polytropic equation of state with adiabatic index 
7 = 7/5. In our spherical code we calculate the shock tube in a spherical shell with thickness d = A 
cm at a large radius of tq = 10000 cm in order to approximate slab symmetry. In the initial state, 
the discontinuity is placed at vq. The left hand side state is {u,p,e,p}j^ = {0,1,2.5,1} and the 
right hand side state is {u, p, e,p}^ = {0, 0.125, 2, 0.1} in cgs units. The discontinuity in the initial 
profile of the variables y is smoothed around tq by a hyperbolic tangent function 

y{r) = y{< ro) + ^ {y{> ro) - y{< ro)) (1 + tanh {g (r - ro))) (74) 

whith a profile slope of 5 = 150 at the edge ro. The comparison of the numerical and the analytical 
results is shown in Fig. (4). As expected from the conservative formulation, the shock speed 
exactly matches the analytical solution and the jump conditions arc fulfilled. Somewhat smoother 
transitions in the numerical calculation around the contact discontinuity and the rarefaction wave 
are on the one hand due to the smoothened initial state and on the other hand to the diffusion 
according to equation (33) that is given by the advection scheme. 

In order to test also the relativistic regime, we set up initial conditions for a relativistic 
shock. The initial state is {u, p,e,p}j^ = {0, 1,2.5 x 10^^,10^^} for r < ro and {u, p,e,p}^ = 
(0,0.125,2 X 10^^,10^^} for r > tq. The velocity u reaches about half of the velocity of light. 
This is much more than expected in any supernova calculation and enough to distinguish rela- 
tivistic effects in the result. In the reproduction of the exact solution, we followed the approach 
of Centrella &: Wilson (1984) and replaced their numerical integration of equation (C13) with an 
analytical solution^. Analytical solutions to the special relativistic shocktube problem in an inertial 
coordinate frame were also derived by Marti Sz. Miiller (1994). Note, however, that the comoving 
coordinates used in our code do not limit to a global inertial frame in the special relativistic limit. 



^To ease the search for an analytical solution to the equivalent of equation (C13) in (Centrella & Wilson 1984), 
we introduce the abbreviations y = TE/D and a:: = (1 — V^)/{V — 5), where T, E, D, V, and ^ have been defined 
as described in (Centrella & Wilson 1984). They are not to be confused with the notation in this paper. Equation 
(C12) can then concisely be written as 

= (. + V)^-^. (75) 
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radius r [cm] radius r [cm] 



Fig. 4. — Tliis figure sliows a siiocktube calculation with initial data leading to nonrelativistic 
velocities. The dots indicate the location of the grid points in the numerical solution. The solid 
line represents the analytical solution. In a time dependent visualization of this graph one would 
observe the grid points moving in accordance with the shock front. It is mainly the space coordinate 
of a grid point in the shock front that changes value, and not so much the represented physical 
state. 
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While the areal radius, r, the enclosed rest mass, a, and the velocity, u, share their definition 
in both, the comoving Misner-Sharp coordinates and the inertial Schwarzschild coordinates, the 
definition of coordinate time is not the same in these coordinate choices. The correspondence is 
straightforward for the region at rest in the inertial frame and the region containing the matter 
between the contact discontinuity and the shock front because it is instantaneously accelerated to 
constant speed as the shock passes through. The fluid elements in the rarefaction wave have a more 
complicated velocity history and show an individual time lapse with respect to the inertial frame. 
In Fig. (5) we compare the numerical solution to the analytical solution. This test demonstrates 
the accurate implementation of the special relativistic conservation laws that lead to the correct 
relativistic shock jump conditions. 



4.3. Oppenheimer-Snyder dust collapse 

The collapse of a uniform and pressureless dust cloud was first calculated in general relativity 
by Oppenheimer & Snyder (1939). Although most of the relativistic terms vanish for negligible 
pressure, this test problem is often recommended for general relativistic hydrocodes. It mainly 
tests the implementation of the momentum equation, because each fluid element behaves like a 
free particle in the self- gravitating pressureless gas. Initial configuration is a dust cloud of two 
solar masses with a density of /? = lO^gcm^'^ In order to assure a negligible pressure, we set the 
initial temperature in the dust cloud to T = lO^^MeV". We had to switch off the adaptive grid for 
this run because spurious variations in the temperature profile lead to grid motions with advection 
inaccuracies that showed a positive feedback on the temperature variations. The result presented 
in Fig. (6) demonstrates that the match of the lapse function to the Schwarzschild metric outside 
of the cloud and that the discretization of the momentum equation are fine. 



4.4. Homologous collapse 

Here we present the homologous collapse that was first analyzed by Goldreich &: Weber (1980). 
Again we use a polytropic equation of state with p = pT/nib and e = (7 — l)~^T/mb, where T is 
the temperature and 7 = 4/3 the adiabatic index. With this equation of state, p = (7 — l)/oe, one 



Instead of dD/dV, leading to equation (C13), we rather evaluate 



After the substitution of z for the variable y according to the relation y = l/2{cosh{z) — 1), we find with dy = 
^y{y+l)dz the ODE 

dV i-v ^ ' 

that has the analytical solution V = tanh ^ ^r_i ^ ■ 
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Fig. 5. — This figure shows a shocktube calculation that leads to relativistic velocities. The physical 
velocity f = u/7 is large enough to distinguish relativistic effects. The dots indicate the location 
of the grid points in the numerical solution. The solid line represents the analytical solution. Note 
that the time measured in our comoving coordinates differs from the time measured in an inertial 
frame. 
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Fig. 6. — This figure shows the collapse of a homologous dust cloud in comoving coordi- 
nates. The numerical solution (dots) can hardly be distinguished from the analytical solu- 
tion (solid lines) at early times. The lines represent mass traces at enclosed baryon masses of 
[0.013, 0.266, 0.519, 0.772, 1.03, 1.28, 1.53, 1.78] solar masses. 
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easily derives the time derivative 



dt \p-y 



) 



7- 1 Tde d /1\ 
p^-i [di ^^di \p) 



= 0. 



(78) 



Thus, the quantity p/p^ is expected to remain constant during the dynamical evolution. Our 
finite difference representation is chosen to conserve total energy to machine precision. The check 
for constant entropy in an adiabatic evolution therefore provides important information about the 
accuracy of the evolution of the internal energy. We run this test with the same artificial viscosity as 
in the realistic simultions in order to check its (non-)influence during the homologous compression 
in the collapse phase. Following Yamada (1997) we start with a hydrostatic isentropic state with 
central density pc = 10^ g cm~^ and temperature Tc = 0.2 MeV. This configuration is uniformly 
cooled until we achieve marginal stability. We then reduce the pressure further by a total of 3% 
(Goldreich k. Weber 1980), distributed over 300 time steps. This leads to a smooth transition 
into collapse without initiating pressure waves that disturb the infall of the surface and lead to 
dissipation of kinetic energy in small shocks. We update the surface pressure proportionally to the 
evolution of the central pressure in order to prevent the surface from drifting out of the homologous 
motion owing to the lack of an outer part of the star. The test is carried out in the Newtonian 
limit, the fully relativistic run shows a less homologous velocity profile, but similar deviations in 
the enetropy. The homologous collapse of an isentropic polytrope is shown in Fig. (7). Snapshots 
are taken for every decade in the central density. The radius of the profiles is normalized to the 
surface radius. Graph (a) shows the entropy evolution. The uppermost solid line belongs to the 
initial state. The lowermost dash-dotted line belongs to a central density of 10^ g cm~^ and clearly 
monitors the 3% pressure reduction applied to initiate the collapse. The entropy deviation grows 
with increasing compactness, but does not exceed the 3% level until nuclear density is achieved in 
the center of the star. Graph (b) shows the nicely homologous velocity profiles with the velocity 
being proportional to the radius. Graph (c) and (d) visualize the selfsimilar evolution of density 
and temperature. 



Finally, we compare the shock tube result calculated with AGILE to the shock tube data pro- 
vided by Raph Hix and Alan Calder from the EVH-1 implementation of PPM. Piecewise Parabolic 
Methods became a standard in hydrodynamics codes while conventional finite difference schemes 
gathered some patina, illustrated by their dependence on the unspeakable artificial viscosity. In 
Fig. (8), a comparison between the two methods is shown and related to the analytical solution. 
The solution of EVH-1 was calculated on 100 equidistant grid points while AGILE used the same 
number of adaptive grid points. The reconstruction routines in EVH-1 do not implement a contact 
discontinuity steepening. A fiattening parameter / < 1 is introduced to provide additional dissipa- 
tion in the shock if a relative pressure jump exceeds e = 0.33 as defined by Colella & Woodward 
(1984) in Eq. (A.l). Following their prescription given in Eq. (A. 2), this flattening parameter has 



4.5. Comparison to PPM 
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Fig. 7. — The homologous collapse of an isentropic polytrope is shown. Snapshots are taken 
when the central density reaches [10^, 10^, 10^°, 10^^, 10^^, 10^^, 10^^] g cm~^. Graph (a) shows the 
entropy evolution. The uppermost solid line belongs to the initial state. The lowermost dash- 
dotted line belongs to a central density of 10^ g cm~^ and clearly shows the 3% pressure reduction 
applied to initiate the collapse. The entropy deviation grows with increasing compactness. Graph 
(b) shows the velocity profiles. Graph (c) and (d) visualize the selfsimilar evolution of density and 
temperature. 
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Fig. 8. — A nonrelativistic shocktube calculation comparing the results of AGILE and EVH-1. 
Both codes use 100 grid points. 



-34- 



been set to 

with w(^) = 0.75 and w^^) = 5. In AGILE, we have set a relatively small artificial viscosity with 
AZ = 2 X 10~^ cm in order to bring the misunderstanding about accuracy and artificial viscosity 
into prominence. The adaptive grid technique renders an unrivalled resolution in the shock front. 
A simulation on an equidistant grid would require 32000 grid points to achieve a similar resolution 
at the shock. Moreover, the shock position and plateaus in the physical state are accurate. With 
a conservative formulation of the equations, the implicit algorithm solves numerically for the same 
jump conditions at the shock front as are solved analytically in a Riemann solver. The artificial 
viscosity only affects the shock width and provides the means for energy dissipation in the shock 
that is always present in a realistic fluid, even if the shock spreads over a narrower spatial range. 
The weakness of the adaptive grid, however, is the rarefaction wave. The diffusivity induced by 
advection on the dynamical adaptive grid as discussed in section 3.2 is responsible for the more 
indistinct rarefaction wave and contact discontinuity. According to equation (33), the artificial 
diffusion is large at increased grid spacing and high grid velocities. As expected, the PPM result 
shows less difficulty in the rarefaction wave, owing to its more uniform grid point distribution 
and smaller diffusivity. The run with EVH-1 required 40 time steps with 4 x 10~^ CPU seconds 
per explicit time step on an IBM RS/6000 SP, the run with AGILE required 429 time steps with 
0.0713 CPU seconds per implicit time step. 295 time steps were used until physical time 0.35 s 
and 134 time steps were needed from 0.35 to 0.7 s. The efficiency increases in the second part 
of the simulation because, as illustrated in Fig. (2), the adaptive grid can take advantage of the 
self-similarity of the flow. About 30% of the computational effort is spent with the inversion of 
the Jacobian. There are many ways to compare the performance depending on the importance of 
specific features. If one would request the same maximum resolution, the number of time steps 
required in EVH-1 would dramatically scale up by a factor of 320 due to the Courant restriction 
on the time step. It is also important to note that the CPU time per time step is not relevant for 
the realistic applications we have in view. The computational effort is dominated by the radiation 
transport or the nuclear reaction network. Important is the overall number of zones and the ability 
to allocate them in regions with interesting physics. Moreover, the total computational effort will 
be proportional to the number of time steps. In the shock tube calculation, the Courant condition 
of explicit schemes is only challenged by resolution. In the supernova application, the restriction 
is much more severe in the compact proto-neutron star, where the sound speed quickly reaches a 
third of the velocity of light. 



5. Conclusion 

We report on AGILE, a first order implicit solver for stiff algebro-differential equations. AGILE 
builds the Jacobian by automatized numerical finite differences. It is therefore flexible with respect 
to changes and extensions in the physical equations and robust against programming errors in 
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the now obsolete implementation of the Jacobian. As discussed in the introduction, implicit finite 
differencing is especially useful in astrophysical applications where the characteristic time scale of 
physical processes varies by orders of magnitude. 

In spherical symmetry, we have written the Einstein equations in comoving coordinates and 
conservative form. These equations are implemented in a module for AGILE to obtain a spherically 
symmetric general relativistic hydrodynamics code. We recognize that a dynamical adaptive grid 
is nothing else than a specific choice of shift vectors in the view of a general relativistic 3+1 
decomposition. The adaptive grid technique of Winkler, Norman, k, Mihalas (1984) and Dorfi & 
Drury (1987), that was developed for Newtonian hydrodynamics, has been extended to relativistic 
space-time and used as a recipe for a stable implementation of shift vectors that focus the numerical 
work to physically important regions in the simulation. 

With an application in stellar core collapse and supernova simulations in view, special care of 
energy conservation was taken. Wc present the detailed wrinkles in our finite difference representa- 
tion of the physical equations. Our guidelines were (i) a consistent discretization of interdependent 
constituents to achieve an exact implementation of physical conservation laws, (ii) avoidance of 
cancellation errors by the independent evolution of small quantities, (iii) the absorption of numeri- 
cal defects into the physical concept of viscosity for an overall consistency. Thus, the discretization 
of the total energy equation is matched to the finite difference representation of the momentum 
and internal energy equation and leads to energy conservation at machine precision. We quantify 
the diffusion introduced by advection in the finite difference representation of the momentum equa- 
tion. Energy conservation requires a compensating pressure term in the energy equation that we 
consistently absorbe into the concept of artificial viscosity. Our artificial viscosity is based on the 
tensor viscosity approach of Tscharnuter &; Winkler (1979) that was extended to general relativity 
by direct inclusion of the viscous stress-energy tensor into the derivation of the Einstein equations 
(Liebendorfer, Mezzacappa, & Thielemann 2001). 

We achieve accurate results in the standard test problems for the supernova application, i.e., 
Sedov point blast explosion, nonrelativistic shock tube, relativistic shock tube, Oppcnheimer-Snyder 
dust collapse, and in the homologous collapse. We also compare our shock tube solution to a 
simulation with the Piecewise Parabolic Method and demonstrate a superior resolution of the 
shock with the adaptive grid. The overall reduction of computational zones, without compromise 
in the resolution of physically important regions, is especially helpful if the hydrodynamics code 
is coupled to computationally more expensive processes, as e.g. Boltzmann radiation transport 
or a nuclear reaction network. The comparison with PPM confirms that the artificial viscosity in 
the adaptive grid technique does not affect the accuracy of the shock speed and the fulfillment of 
the correct jump conditions. However, with the low order advection scheme chosen in this work, 
considerable diffusion happens in the rarefaction wave where the adaptive grid does not opt for a 
high resolution. 

A previous version of AGILE has been used to simulate the hydrodynamical part in X-ray 
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burst simulations on an accreting neutron star (Rembges 1999). An AGILE module implementing 
general relativistic hydrodynamics in polar slicing and radial gauge has been tested and extended 
by Miiller (2000) to study the effect of dark matter in the formation of supermassive black holes 
in the early universe. A Newtonian module in Lagrangian coordinates has been used in stellar core 
collapse and postbounce evolution simulations with complete 0(v/c) Boltzmann neutrino transport 
(Mezzacappa et al. 2001). The module described in this paper, in combination with general 
relativistic Boltzmann neutrino transport, lead to the general relativistic simulation of stellar core 
collapse and postbounce evolution (Liebendorfer et al. 2001). 
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